Differentially Expressed Genes and Molecular Susceptibility to Human Age-Related Diseases

Mainstream transcriptome profiling of susceptibility versus resistance to age-related diseases (ARDs) is focused on differentially expressed genes (DEGs) specific to gender, age, and pathogeneses. This approach fits in well with predictive, preventive, personalized, participatory medicine and helps understand how, why, when, and what ARDs one can develop depending on their genetic background. Within this mainstream paradigm, we wanted to find out whether the known ARD-linked DEGs available in PubMed can reveal a molecular marker that will serve the purpose in anyone’s any tissue at any time. We sequenced the periaqueductal gray (PAG) transcriptome of tame versus aggressive rats, identified rat-behavior-related DEGs, and compared them with their known homologous animal ARD-linked DEGs. This analysis yielded statistically significant correlations between behavior-related and ARD-susceptibility-related fold changes (log2 values) in the expression of these DEG homologs. We found principal components, PC1 and PC2, corresponding to the half-sum and the half-difference of these log2 values, respectively. With the DEGs linked to ARD susceptibility and ARD resistance in humans used as controls, we verified these principal components. This yielded only one statistically significant common molecular marker for ARDs: an excess of Fcγ receptor IIb suppressing immune cell hyperactivation.

This could be why mainstream transcriptome-profiling studies of ARD susceptibility versus ARD resistance in human volunteers [33,[36][37][38][39][40][41][42][43] and animals [30,32,35, are focused on the differentially expressed genes (DEGs) that are specific to gender, age, tissues, and pathogeneses. All previous research contributes to progress in predictive, preventive, personalized, participatory (4P) medicine [65] and helps us understand where, how, why, when, and what disorders can affect people depending on their genetic background, medical history, and lifestyle. Because none of us is able to escape an ARD, we expected that a meta-analysis of all available ARD-linked DEGs would eventually reveal the most common, universal theranostic molecular marker that will be permanently available in any tissue of anyone's organism.
Because hypertension contributes to vascular aging [66] and vice versa, we have previously studied, within this mainstream paradigm, the inbred ISIAH (Inherited Stress-Induced Arterial Hypertension) rat strain [67] and sequenced transcriptomes in the brain stem [47], hypothalamus [48], renal medulla [49], renal cortex [50], and adrenal glands [51], with WAG rats used as controls. Additionally, we have recently obtained and compared transcriptomes of the midbrain tegmentum [68] in gray rats of a tame and an aggressive outbred strain [69][70][71][72]; this allowed us to identify, in in silico settings, potential molecular markers for neoteny, which has a power to reverse ARDs [73]. Furthermore, we have recently profiled transcriptomes in the hippocampus [74] of tame versus aggressive rat strains, in which deficient β-protocadherins and β-hemoglobin were found to be the statistically significantly most common theranostic molecular markers for ARD-related hypertension. That is why we sequenced, in the same way, one more transcriptome, that of the periaqueductal gray matter (PAG), in the tame versus aggressive rats, identified the corresponding behavior-related PAG-associated DEGs in them, and compared these DEGs with their homologous PubMed-based [75] ARD-linked DEGs in animals [30,32,35,68,74] and humans [33,36-43] to find out whether there are invariant molecular markers for such diseases among them.

RNA-Seq and Mapping to the Reference Rat Genome
We focused on the PAG because the activity of this brain structure contributes to elevated pain tolerance [76], sociability, intelligence (IQ), and humor processing with increasing age [77]. In fact, our best hope was that if we had had all possible ARD molecular markers at hand, we would have eventually understand which of them might be the best helpers to relieve the suffering of ARD patients. We profiled the PAG transcriptome of three tame adult male gray rats (Rattus norvegicus) against three aggressive conspecifics on an Illumina NextSeq 550 system (see Section 4.2). The rats came from two outbred strains, one tame and one aggressive, maintained at the Institute of Cytology and Genetics of the Siberian Branch of the Russian Academy of Science [72,78] for more than 90 generations using the glove test [79]. The rats were not consanguineous (see Section 4.1). First, we sequenced 210,128,758 reads each 75 nt in length and deposited them in the NCBI SRA database [80] (ID PRJNA668014) (see Table 1). Next, we chose from among them 177,608,837 raw reads (84.5%) via mapping to the reference rat genome RGSC Rnor_6.0, UCSC Rn 6 July 2014 (Table 1). Then we identified 14,039 genes expressed in the PAG of the studied rats. Finally, we selected 39 DEGs using Fisher's Z-test with the Benjamini correction for multiple comparisons and discarded hypothetical, tentative, predicted, uncharacterized, or non-protein-coding genes to minimize the false-positive error rates (Tables 1 and 2). Note: Hereinafter, log 2 : the log 2 -transformed fold change (i.e., the ratio of the expression level of a given gene in tame rats to that in aggressive rats); p and P ADJ : statistical significance according to Fisher's Z-test and the Benjamini correction for multiple comparisons.

Quantitative PCR (qPCR)-Based Selective Verification of the Novel PAG-Related DEGs of the Tame and Aggressive Rats
From the same two strains of rats, we took 16 additional unrelated animals: eight tame and eight aggressive, each scoring 3.5 and -3.5, respectively, on a scale spanning from -4 (most aggressive) to +4 (least aggressive) in the glove test [79] run 1 month before the PAG specimens were sampled (Table 3). Next, from among the 39 DEGs shown in Table 2, we selected Ascl3 and Defb17 (for details, see Table 3). Table 3 contains our qPCR data on these genes in the PAG of the tame and aggressive rats (see Section 4.4). These PCR data appear as the arithmetic mean ± standard error of the mean of Ascl3 and Defb17 expression levels normalized to those of three reference genes, B2m (β-2-microglobulin) [81], Hprt1 (hypoxanthine phosphoribosyltransferase 1) [82], and Rpl30 (ribosomal protein L30 [83]), in triplicate, according to the guidelines [84]. The rightmost column of Table 3 shows arithmetic mean estimates of the expression levels of Ascl3 and Defb17 in the PAG of the tame and aggressive rats.
As can be seen from Figure 1a, both Ascl3 and Defb17 are overexpressed in the PAG of the tame rats (white bars) compared to the aggressive rats (gray bars) according to our qPCR data, this overexpression being significant (p < 0.01, double asterisk) according to both the nonparametric Mann-Whitney U test and the parametric Fisher's Z-test. character "**") denotes statistical significance at p < 0.01 according to nonparametric Mann-Whitney U test and parametric Fisher's Z-test. (b) Pearson's linear correlation between the relative expression levels of Ascl3, Defb17, and the reference genes [B2m (β-2-microglobulin) and Rpl30 (ribosomal protein L30) appearing as circles on the plot] in the tame versus aggressive rats is statistically significant, whether measured experimentally using RNA-Seq (x-axis) or qPCR (y-axis) and expressed on the log 2 scale (see "Materials and Methods"). Dashed and dot-dash lines denote linear regression and its 95% confidence interval boundaries calculated using STATISTICA (Statsoft TM , Tulsa, OK, USA); r and p are Pearson's linear correlation coefficient and its statistical significance, respectively. This finding independently supports our RNA-Seq data ( Table 2). Finally, Figure 1b shows Pearson's linear correlation (r = 0.997, p < 0.005) between the log 2 values ("log 2 " hereinafter means "the log 2 -transformed ratio of the expression level of a given gene between tame and aggressive rats under the given experimental conditions") for four genes, Ascl3, Defb17, B2m, and Rpl30 (all appearing as open circles), within the RNA-Seq (x-axis) and qPCR (y-axis) datasets, both obtained here independently. This is one more piece of independent evidence in support of the relevance of our RNA-Seq data ( Table 2).

Comparison of Known Animal ARD-Linked DEGs with Their Homologs among 39 Novel PAG-Related DEGs of the Tame and Aggressive Rats
We retrieved all transcriptomes of animals with ARD susceptibility and ARD resistance from the PubMed database [75] and collected 43 animal-based human ARD models (Table 4).   Figure 2 shows our procedure of comparison of 37,834 animal ARD-linked DEGs (Table 4) with the 39 PAG-related DEGs of the tame and aggressive rats ( Table 2). In this figure, a Venn diagram and a table explain the procedure that led us to the following conclusion: there is a correlation between the expression of homologous genes in the tame rats and the ARD-susceptible animals. At the first step, we compiled 459 pairs of homologous DEGs, where one of 39 DEGs was taken from Table 2 and its homolog was chosen from among 37,834 DEGs in Table 4, both now being present in Table S1 (see "Supplementary Materials"). At the second step, we for the first time found the following statistically significant correlations between behavior-related and disease-susceptibility-related log 2 values of the homologous animal DEGs: Pearson's linear correlation (r = 0.13, p < 0.01), Spearman's rank correlation (R = 0.14, p < 0.005), Kendall's rank correlation (τ = 0.10, p < 0.0025), and the Goodman-Kruskal generalized correlation (γ = 0.10, p < 0.0025) ( Figure 2). These correlations suggest a certain similarity in the expression patterns of the novel PAG-related DEGs between the tame and aggressive rats and a similarity in the expression patterns of their homologous DEGs between ARDsusceptible and ARD-resistant animals. The biological sense of these similarities was elucidated at the final step. We processed entries in Table S1 with principal component analysis in the bootstrap mode with the freely available toolbox PAST4.04 [85]. In the result, we found two principal components, major PC1 and minor PC2, corresponding to the half-sum and the half-difference of the behavior-related and ARD-susceptibility-related log 2 values for the homologous animal DEGs, respectively. Here, the half-sum of the behavior-and ARD-associated log 2 values, which appears as the major PC1, implies that the sense of the significant positive correlations found between these log 2 values at step 2 is that their signs in the novel PAG-related rat DEGs and in the known ARD-linked animal DEGs are matching. Their formulas and 95% confidence intervals are given in the bottom part of Figure 2.  Table 4; the xand y-axes correspond to columns iii and x of Table S1; circles correspond to rows of Table S1; R, τ, and γ are the coefficients of Spearman's rank correlation, Kendall's rank correlation, and the Goodman-Kruskal generalized correlation, respectively, calculated using STATISTICA (Statsoft TM , Tulsa, OK, USA); PC1 and PC2: principal components calculated in the bootstrap-based refinement mode of the PAST4.04 software [85].

Animal ARD-Linked DEGs Are Relevant to Humans
We additionally searched for all the PubMed DEGs in ARD patients and otherwise healthy volunteers [75] (see Table 5). The total number of human ARD-linked DEGs found is 14,535 in ten tissues from 14 binary "susceptibility versus resistance" models of human ARDs (the rightmost column of Table 5) [33,36-43]. Figure 3 illustrates a modification of the previous procedure, with the animal ARD-linked DEGs replaced by their human counterparts (Table 5). Table 5. DEGs in the binary "susceptibility versus resistance" models of human ARDs (PubMed data, [75] Table 2. The results serve as independent control medical data (Table S2). Although no correlation was found between the behavior-related log 2 values for the novel PAG-related rat DEGs and the ARD-susceptibility-related log 2 values for their human homologs, their half-sum and half-difference corresponded to principal components PC1 (major) and PC2 (minor) within their quite similar 95% confidence intervals (Figure 3). This finding confirmed the results of our meta-analysis of ARD susceptibility versus ARD resistance in animals.

Searching for ARD Molecular Markers among Human Genes Orthologous to 39 Novel PAG-Related DEGs of the Tame and Aggressive Rats
First, we used the PubMed database [75] and characterized each of the 39 novel PAGrelated DEGs of the tame and aggressive rats ( Table 2) by answering the question as to how under-or overexpression of their human orthologs can aggravate or alleviate ARDs   (Table S3). Next, for each novel PAG-related rat DEG (Table 2), we determined the sign of the behavior-related log 2 value and found how many of their homologs in the ARD-susceptible and ARD-resistant subjects have matching signs of the ARD-susceptibility-related log 2 value (n = N PC1 ) and how many have opposite signs (n = N PC2 ). This is because the principal components PC1 and PC2 correspond to the half-sum and the half-difference, respectively, of these log 2 values (Figures 2 and 3). Table 6 shows these figures (N PC1 and N PC2 ) alongside their statistical significance estimates according to the binomial distribution both before (p-values) and after (P ADJ values) Bonferroni's correction for multiple comparisons. Table 6. Searching for ARD molecular genetic markers among human genes homologous to 39 novel PAG-related DEGs of the tame and aggressive rats. The number of these homologous DEGs linked to ARD susceptibility and ARD resistance in humans and animals is taken into account. As can be seen from this table, only one of the 39 novel DEGs in the PAG of the tame and aggressive rats is linked with PC1, namely Fcgr2b encoding Fcγ receptor IIb, which suppresses the hyperactivation of immune cells [199] (Table 7). Entries in this table suggest that the immunoregulatory DEGs homologous to rat Fcgr2b are significantly overexpressed in ARD-susceptible human and animal subjects compared to their ARD-resistant peers [35,39,43,51,55,56,58]. This leads us to propose that an excess of the immunoregulators homologous to rat Fcgr2b might be regarded as a candidate theranostic molecular marker for human ARDs. The following facts and hypotheses were thought to be relevant to this question: (i) the number of human infections increases when new wild animals are domesticated and become a bridge between the anthropogenic environment and wildlife [200], (ii) stray cats live longer than pet cats due to stronger disease resistance [201], and (iii) human evolutionary origin involves self-domestication [202], although this hypothesis is still debatable [203].

Rat Gene
With this information in mind, we searched PubMed [75] and retrieved all available transcriptomes of domestic animals compared to their wild conspecifics ( Table 8). As can be seen from the bottom line of Table 8, RNA-Seq data represent 2866 DEGs in eight tissues of seven domestic animal species compared to their seven wild conspecifics studied in ten original experimental works [68,74,[204][205][206][207][208][209][210][211]. By looking at the entries in the left half of Table 9, it becomes clear why underexpression of the human FCGR2B gene, which is orthologous to the rat Fcgr2b gene (this being the only gene we have identified as a theranostic molecular marker for ARD; Table 9), alleviates such diseases [117,118], while its overexpression aggravates them [119,120]. Table 9. Comparison of the effects of unidirectional changes in the expression (a) of the human FCGR2B gene on ARD severity in humans and (b) of its animal homologs on the microevolutionary events leading to domestic and wild animals.  From among the 39 novel PAG-related rat DEGs (Table 2) and these 2866 known animal DEGs (Table 8), we selected seven animal genes homologous to the rat Fcgr2b gene (the right half of Table 9). Furthermore, we transformed the log 2 value characterizing the animal DEGs into either underexpression or overexpression of these animal genes after the domestic and wild animals split from their most recent common ancestor, according to the commonly accepted phylogeny concept [212][213][214][215][216]) ( Table 9).
The left half of Table 10 summarizes the data shown in Table 9 in the form of a 2 × 2 contingency table, where this difference is statistically significant according to both Pearson's χ 2 test (p < 0.01) and Fisher's exact test (p < 0.05). Table 10. Correlations between the effects of unidirectional changes in the expression (a) of the human FCGR2B gene on ARD severity in humans and (b) of its animal homologs on the microevolutionary events leading to domestic and wild animals.

Pearson's χ 2 Test Fisher's Exact Test Alleviating (←) Aggravating (→) χ 2 p
Effect of changes in the expression of animal homologs to human FCGR2B wild 6 1 7.14 10 −2 0.05 domestic 1 6 In this regard, upregulation of the immunoregulators homologous to the rat Fcgr2b gene can serve as a theranostic ARD molecular marker permanently available in any tissue of anyone's organism.

Overexpression of Human Genes Homologous to the Rat Fcgr2b Gene Identified as A Molecular Marker for ARDs Is Consistent with Overexpression of Known ARD-Linked DEGs
We have found that upregulation of the human immunoregulators homologous to the rat Fcgr2b gene can serve as a theranostic ARD molecular marker permanently available in any tissue of anyone's organism (see Table 7). As can be seen from Table 7, the immunoregulators homologous to the rat Fcgr2b gene are excessive in all tissues of all ARD-affected subjects (humans and animals) with only one exception: FCGR1B deficiency in the lungs of a patient with fibrosis in pulmonary arterial hypertension [39] ( Table 7: line #12). This implies that an excess of the human immunoregulators homologous to the rat Fcgr2b gene identified in this work as being a molecular marker for age-related human diseases is consistent with the ARD-linked DEGs found independently by other authors.

Excess of Rat Fcgr2b Identified as A Molecular Marker for ARDs Is Consistent with Independent RNA-Seq Data on the Rise in the Levels of Murine Fcgr2b in Astrocytes with Age
We will discuss this result in comparison with independent RNA-Seq data [217] from a human ARD model using mice aged 7 days, 32 days, 10 weeks, 9.5 months, and 2 years. Although murine Fcgr2b was not identified as a DEG within the model in question [217], we have found a statistically significant increase in the expression of this gene in astrocytes (see Figure 4) from the hippocampus (open circles), striatum (gray circles), and cortex (black circles) of mice with increasing age. This could be regarded as an in vivo piece of independent direct evidence in support to our findings in this work.

Stabilizing Selection Preserves the Expression Pattern of the Human FCGR2B Gene Orthologous to the Rat Fcgr2b Gene, a Molecular Marker for ARDs
Because it now seems interesting to discuss how such an extraordinary expression pattern-that is, with levels increasing with age-can have been preserved, we will examine the proximal promoter of the human FCGR2B gene (see Figure S1). As can be seen from Figure S1a, the current build (No. 155) of the dbSNP database [218] contains only one single-nucleotide polymorphism, SNP rs780467580, within the 70 bp proximal promoter of the human FCGR2B gene. We have previously used our publicly available development, SNP_TATA_Comparator [219], and manually analyzed 15,080 SNPs within 1585 proximal promoters, each 70 bp in length, located upstream of the transcription start sites of proteincoding transcripts from 534 human genes (for review, see [220]). The number of SNPs within these promoters varied from one to 64 (9.51 ± 0.21; mean ± SEM), with the 95% confidence interval being between 9 and 10; therefore, the existence of only one SNP rs34166473 within the examined region of the human FCGR2B gene promoter seems to suggest a statistically significant loss of SNPs (p < 0.01, binomial distribution). It is generally recognized that the biological function of any genomic region containing such a small number of SNPs is of vital importance, and so it is preserved by natural selection [221]. This observation is consistent with the neutrality of this single SNP rs34166473 within the human FCGR2B gene promoter in question, the neutrality having been confirmed with the use of our previously developed publicly available web service SNP_TATA_Comparator [219] (see Figure S1).

The Use of Tame and Aggressive Rats as Belyaev's Laboratory-Animal-Based Domestication Model can Represent an Adequate Domesticated-Animal-Based Model of Human ARDs
In searching for ARD molecular markers, we used tame and aggressive rats as Belyaev's laboratory-animal-based domestication model [69][70][71][72], because we have already found molecular markers for hypertension (incidentally, this being an ARD) and potential molecular markers for neoteny (this can reverse ARDs [73]), which correspond to the hippocampal [74] and midbrain tegmental [68] transcriptomes of these rat strains. As can be seen from Tables 9 and 10, within tissues of the domesticated compared to wild animals there is a statistically significant excess of immunoregulators homologous to rat Fcgr2b according to two independent criteria, Pearson's χ 2 test (χ 2 = 7.14; p < 0.01) and Fisher's exact test (p < 0.05). This implies that domesticated animals are susceptible to ARDs, while their wild conspecifics are not. This result is consistent with the findings reported by Fallahshahroudi and co-workers [201] that, in the wild, natural selection on animals seeks to eliminate affected individuals, while artificial animal selection during domestication does not, as it serves human needs (for example, broilers are the cheapest poultry meat, no matter whether they are susceptible or resistant to ARDs [30,34]). This is in line with our whole-genome observations suggesting that the TATA-binding protein binding sites within the gene promoters in domesticated animals have significantly more candidate SNP markers for rheumatoid polyarthritis [222] and reproductive disorders [223,224], both being ARDs, compared to their wild conspecifics. Overall, we can conclude that Belyaev's model of domestication using tame and aggressive rats as laboratory animals seems to be applicable as one of the adequate animal-based human ARD models.

Our Focus on the PAG of Tame and Aggressive Rats Is Adequate for the Search for the ARD-Linked Rat DEGs
In this work, we focused on the PAG of the tame and aggressive rats because the activity of this brain structure contributes to elevated pain tolerance with age [77]. Our interest was to find ARD-linked rat DEGs such that, if their human homologs were taken into account, we could help reduce the suffering of patients with such diseases (see Section 2.1). Recently, a meta-analysis of the microarray datasets GSE24982, GSE63442, and GSE63651 (from the Gene Expression Omnibus (GEO) database [225]) identified the rat Fcgr2b gene as one of the six most likely hub genes responsible for neuropathic pain and aging [226], just as we found the same rat Fcgr2b gene to be a molecular marker for ARDs (see Table 2, Table 6,  Table 7, Table 9, and Table 10). Apparently, this independent microarray meta-analysis result [226] favors the adequacy of our focus on the PAG of the tame versus aggressive rats in the search for ARD-linked rat DEGs.
Briefly, we used the PAG of the tame and aggressive rats and found rat Fcgr2b overexpression to be a molecular marker for elevated neuropathic pain tolerance in ARDs.

Human Immunoregulatory Genes Homologous to the Rat Fcgr2b Gene Identified as a Molecular Marker for Pain Tolerance in Age-Related Diseases Are Young on the Molecular-Evolution Scale
Because the human FCGR2B gene orthologous to the rat Fcgr2b gene identified as a molecular marker for elevated neuropathic pain tolerance in age-related diseases seems to be under stabilizing selection ( Figure S1), we used our Orthoscape plug-in [229,230] within the Cytoscape software suite [231] and estimated BLAST-based [232] phylostratigraphic age indexes (PAIs) for (a) 39 human genes homologous to the behavior-related PAG-associated rat DEGs identified in this work ( Table 2) and (b) 48 human FCGR2B-associated genes ( Figure 5). The results are in the top and bottom parts of Table S4. As can be seen from the top part of this table, FCGR2B is one of the two youngest behavior-related human genes. Furthermore, the top ten youngest genes among all 48 human genes in its bottom part contain all six human FCGR1A, FCGR2A, FCGR2B, FCGR2C, FCGR3A, and FCGR3B genes, which is a statistically significant event according to the binomial distribution (p < 0.0001). Finally, as can be seen from Figure 6, the 48 FCGR2B-associated human genes examined are significantly younger than the 39 behavior-related human genes, according to the nonparametric Mann-Whitney U test (p < 0.01) and the parametric Fisher's Z-test (p < 0.01). Figure 6. A comparison between 48 human FCGR2B-associated genes using the publicly available toolbox ANDSystem [227] (Figure 5) and the human genes homologous to 39 novel PAG-related rat DEGs (Table 2). Legend: the double asterisk (i.e., double character "**") denotes statistical significance at p < 0.01 according to nonparametric Mann-Whitney U test and parametric Fisher's Z-test; PAI, a gene's phylostratigraphic age index evaluated against the BLAST-based scale [232] using the freely available web service Orthoscape [229]; BLAST-based PAI scale: 0, Cellular organisms; 1, Eukaryota; 2, Opisthokonta; 3, Metazoa; 4, Eumetazoa; 5, Bilateria; 6, Deuterostomia; 7, Chordata; 8, Craniata; 9, Vertebrata; 10, Gnathostomata; 11 Human immunoregulatory genes homologous to the rat Fcgr2b gene that we have identified as a molecular marker for elevated neuropathic pain tolerance in ARDs are expressed at increased rates mostly in conventional, monocyte-derived, and plasmacytoid dendritic cells as well as macrophages [233], which occur in most tissues, where they are critical to tissue homeostasis [234]. According to the GeneCard database [235], a large number of microarray, RNA-Seq, and proteomics experiments have detected molecular products expressed from the human genes FCGR1A, FCGR2A, FCGR2B, FCGR2C, FCGR3A, and FCGR3B in the majority of the human tissues studied. Because disruptions in the cellular-senescence-associated tissue homeostasis compromise the correct activation of immune responses to pathogens and cancer cells [236,237], these human immunoregulatory genes can be permanently overexpressed in any tissue of anyone's organism as molecular markers for elevated neuropathic pain tolerance in ARDs (see Table 7).

Animals
The animals used were adult male gray rats (R. norvegicus) artificially bred for over 90 generations for either aggressive or tame behavior towards humans as two outbred strains. The rats were kept under standard conditions of the Conventional Animal Facility at the ICG SB RAS (Novosibirsk, Russia), as described elsewhere [72,79,238]. The total number of rats was 22 (11 tame and

RNA-Seq
Total RNA was isolated from~100 mg of the PAG tissue samples of tame (n = 3) and aggressive (n = 3) rats using the TRIzol™ reagent (Invitrogen, Carlsbad, CA, USA). The quality of the total-RNA samples was measured on a Bioanalyzer 2100 (Agilent, Santa-Clara, CA, USA). Samples with optimal RNA integrity numbers (RINs) were selected for further analysis. Furthermore, the total RNA was analyzed quantitatively on an Invitrogen Qubit™

Mapping of RNA Sequences to the R. norvegicus Reference Genome
The primary raw Fastq files were examined using a quality-control tool for high-throughput sequencing data (FastQC; https://www.bioinformatics.babraham.ac.uk/projects/fastqc; accessed on 19 December 2018). Next, we improved the quality of the raw reads using the Trimmomatic tool [240] in a step-by-step manner as follows: (i) discarded a base from either the start or end position if the quality was low; (ii) trimmed bases with a sliding-window method, and (iii) eliminated any remaining reads that were less than 36 bases in length. After that, we aligned the trimmed reads to the R. norvegicus reference genome (RGSC Rnor_6.0, UCSC version Rn6, July 2014 assembly) using the TopHat2 toolbox [241]. Next, we reformatted these alignments into sorted BAM files in SAMtools version 1.4 [242]. Then we assigned the reads in question to these genes using the htseq-count tool from the preprocessing software HTSeq v.0.7.2 [243] along with gtf files containing the coordinates of the rat genes according to Rnor_6.0 and an indexed SAM file. Finally, we used DESeq2 [244] via the web service IRIS (http://bmbl.sdstate.edu/IRIS/; accessed on 16 January 2020), rated differential expression levels of the rat genes, and, to minimize false-positive error rates, applied Fisher's Z-test [245] with the Benjamini correction for multiple comparisons as well as discarded all hypothetical, tentative, predicted, uncharacterized, and non-protein-coding genes.

qPCR
To examine independently and selectively the novel tame versus aggressive rat PAG DEGs identified here (Table 2), we performed a qPCR control assay on the total RNA taken from the remaining samples of the PAG of tame (n = 8) and aggressive (n = 8) rats. First, we isolated total RNA using TRIzol™, purified it on Agencourt RNAClean XP Kit magnetic beads (Beckman, #A63987), and quantified it using a Qubit™ 2.0 fluorometer (Invitrogen/Life Technologies) and a Qubit™ RNA High-Sensitivity Assay Kit (Invitrogen, cat. # Q32852). Next, we synthesized cDNA using the Reverse Transcription Kit (Syntol, #OT-1). Finally, we designed oligonucleotide primers for qPCR using the web service PrimerBLAST [246] (Table 11). Table 11. qPCR primers selected via the publicly available web service PrimerBLAST [246].

Human Genes
We analyzed the 39 human genes that are orthologous to the 39 PAG-related DEGs of the tame and aggressive rats ( Table 2). With reliance on PubMed [75], we gathered clinical information about whether downregulation or upregulation of each of these 39 human genes can alleviate or aggravate ARDs (Table 9 and Table S3).

DNA Sequences
For in silico analysis of the human genes encoding the novel candidate molecular markers for ARDs, we retrieved the DNA sequences and SNPs of their 70 bp proximal promoters from the Ensembl database [247] and the dbSNP database [218], respectively, using the UCSC Genome Browser (reference human genome assembly GRCh38/hg38) [248] in both dialog mode and automated mode using the Bioperl toolkit [249] ( Figure S1).

A Knowledge Base for Domestic Animal DEGs Whose Human Orthologs Can Affect ARD Severity
In files with the flat Excel-compatible textual format, we have documented all proposed associations between the domestic and wild animal DEGs homologous to the 39 novel DEGs that we found in the PAG of the tame and aggressive rats. We have also documented how downregulation or upregulation of the human genes homologous to these PAGrelated rat DEGs can affect ARD severity. Finally, we have added our current findings to our previously published public PetDEGsDB knowledge base, its new build being freely available at www.sysbio.ru/domestic-wild (accessed on 14 December 2022) in the MariaDB 10.2.12 database management system (MariaDB Corp AB, Espoo, Finland).

Data Mining of Literature Sources and Databases Publicly Available on the Internet
We conducted data mining using the automated mode of our previously published freely available web service ANDSystem [227], with "Human, FCGR2B, gene, protein" as input keywords, with all the other parameters set at their default values.

In Silico Estimation of the BLAST-Based PAIs of a Given Human Gene
We estimated the BLAST-based [232] PAIs for a given human gene whose NCBI Entrez gene number served as input for our Orthoscape plug-in [229,230] within the Cytoscape software suite [231]. The output was the most recent common ancestor of all the animal species whose DNA sequence of this gene is already known. The following evolutionary rank scale was used: 0,

Statistical Analysis
We performed the Mann-Whitney U test, Fisher's Z-test, Pearson's linear correlation test, Goodman-Kruskal generalized correlation test, Spearman's and Kendall's rank correlation tests, Pearson's χ 2 test, Fisher's exact test, and the binomial-distribution analysis using appropriate options in the standard software STATISTICA (Statsoft TM ). Furthermore, using the PAST4.04 software package [85], we conducted a principal component analysis in the bootstrap refinement mode via its mode selection path "Multivariate"→"Ordination"→ "Principal Components (PCA)"→"Correlation"→"Bootstrap."

Conclusions
First, we have profiled the PAG transcriptome in three tame adult male rats versus three aggressive conspecifics, all animals being unrelated, and made the primary raw reads publicly available (NCBI SRA database ID: PRJNA668014) [80]. With the use of this transcriptome, we have found 39 PAG-related DEGs whose statistical significance (P ADJ < 0.05, Fisher's Z-test with the Benjamini correction for multiple comparisons) was acceptable ( Table 2). We have selectively verified these DEGs with independent experimental analyses (qPCR) of eight other tame and eight other aggressive adult male rats from unrelated litters of the same two outbred strains (Table 3 and Figure 1).
Secondly, we have found that Fcγ receptor IIb overexpression is a statistically significant molecular marker for ARDs. To come to this conclusion, we compared 39 novel DEGs in the PAG of tame and aggressive rats to their known homologs associated with ARDs in animals and humans, using correlation and principal component analysis, as well as Bonferroni's correction for multiple comparisons.